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Abstract. We consider a methodology based in B-splines scaling functions to numerically 
invert Fourier or Laplace transforms of functions in the space L 2 (M). The original function 
is approximated by a finite combination of j th order B-splines basis functions and we provide 
analytical expressions for the recovered coefficients. The methodology is particularly well 
suited when the original function or its derivatives present peaks or jumps due to discontinu- 

' ities in the domain. We will show in the numerical experiments the robustness and accuracy 

. of the method. 
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Fourier and Laplace transforms are useful in a wide number of applications in science and 
engineering. As it is well known, the Fourier transform is closely related to the Laplace 
transform for zero value funcions on the negative time axis. There is a strong interest in the 
efficient numerical inversion of Laplace transforms ([1],[2]) and Fourier transforms, due to the 
fact that the solutions to some problems are known in the transform domain rather than in 
the original domain. An increasing number of papers have recently appeared to invert Fourier 
transforms with wavelets, like [S] and [5] with coiflets wavelets and [TU] with Mexican, Morlet, 
Poisson and Battle-Lemarie wavelets. In particular in the Financial Engineering context, [TT] 
inverts a Laplace transform by means of B-splines wavelets of order 1. 

Recently, a new method called the COS method developed in [7j for solving the inverse 
Fourier integral is capable to accurately recover a function from its Fourier transform in a 
short CPU time being as well of very easy implementation. However, when the function to 
be recovered presents discontinuities or it is highly peaked, a lot of terms in the expansion 
must be considered to reduce the approximation error. 

In this paper we present a novel approximation based on B-splines scaling functions to 
numerically invert Fourier transforms that it is particularly suited for functions that exhibit 
peaks and/or jumps in its domain. There are some properties that, at first glance, make these 
basis functions particularly suited to approximate such non-smooth functions. B-splines are 
^ ! the most regular scaling functions with the shortest support for a given polynomial degree. 

Another important fact is the explicit formulation in the time (or space) domain as well as in 
the frequency domain. However a slight drawback of these basis functions is that the system 
that they form is semi-orthogonal, i.e., the scaling functions are orthogonal among different 
scales but not necessarily at the same scale. 

In previous work [T3] the authors numerically inverted the Laplace transform of a distri- 
bution function in the interval [0, 1] making use of the Haar scaling functions. Now, in the 
present paper, we consider the problem of inverting the Fourier transform to recover an L 2 (IR) 
function by approximating it by a finite sum of B-splines scaling functions of order j (where 
j = is the particular case of the Haar system). We also provide a list of the different errors 
accumulating within the numerical procedure. So, here, the Fourier inversion is carried out 
with B-splines scaling functions rather than with wavelet functions ([10]) or a combination 
of both ([8], [9]). We fix the scale parameter in the wavelet expansion and only remaining is 
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the translation parameter, facilitating the inversion procedure. Furthermore, we provide an 
analytical expression for the coefficients of the approximation. 

As will be shown in the section devoted to numerical examples, the Wavelet Approximation 
(WA) that we present is well capable to detect jumps or peaks produced by discontinuities 
in the function itself or in first derivatives. On the contrary, the COS method is better to 
approximate analytical functions. 

The paper is organized as follows. In Section 1.1 we give a short literature review re- 
garding the Laplace transform inversion. Section 2 gives a brief introduction concerning 
multiresolution analysis and B-splines scaling functions. In Section 3 we present the Wavelet 
Approximation method to recover functions from its Fourier transform by means of B-splines 
scaling functions. Section 4 is devoted to the COS method, while numerical experiments, 
comparing the Wavelet Approximation and the COS method are shown in Section 5. Finally, 
Section 6 concludes. 

1.1. Laplace transform inversion. Suppose that / is a real- or complex- valued function 
of the variable x > and s is a real or complex parameter. We define the Laplace transform 
of / as, 

(1) /(a) = / e- sx f(x)alx = lim / e~ sx f(x)dx, 

whenever the limit exists (as a finite number). 

We state the inverse transform as a theorem (see [6 J for a detailed proof). 

Theorem 1. (Bromwich inversion integral) If the Laplace transform of f(x) exists, then, 

-i ra-\-ik 

(2) /(x)= lim — / f(s)e sx ds, x > 0, 

fc^+oo 2m J a _ ik 

where \f{x)\ < e Sx for some positive real number £ and a is any other real number such that 
a > E. 

The usual way of evaluating this integral is via the residues method taking a closed contour, 
often called the Bromwich contour. 

In this section we present some numerical algorithms to invert the Laplace transform. A 
natural starting point for the numerical inversion of Laplace transforms is the Bromwich 
inversion integral stated in Theorem [TJ If we choose a specific contour and perform the 
change of variables s = cr + m in (j2j), we obtain an integral of a real valued function of a 
real variable. Then, after algebraic manipulation and after applying the Trapezoidal Rule we 
obtain, 

(3) f(x) ~ f h (x) = 3ft (/(a)) + -J— » (7(o- + ikh)) cos(khx), 

k=l 

where 9ft(z) denotes the real part of z. A detailed analysis of the errors can be found in [2]. 

As pointed out in [2], the Bromwich inversion integral is not the only inversion formula 
and there are quite different numerical inversion algorithms. We refer the readers to the 
Laguerre series representation given in [T], which is known to be an efficient method for 
smooth functions. However, if / is not smooth at just one value of x, then the Laguerre 
method has difficulties at any value of x. 

In the context of numerical Laplace inversion, [8] recovers the function / with a procedure 
based on wavelets. They consider s = f3 + iu in expression ([2]), where a> is a real variable 
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and f3 is a real constant that fulfills f(x)e~ l3x G L 2 (R), assuming that f(x) = when x < 0. 
Then, equation ([1]) can be rewritten as, 

/+oo 
e- px e~ iu]X f{x)dx. 
-oo 

Defining, 

(4) h{x) = f{x)e~P x , then h{u) = J(/3 + iu), 

where h denotes the Fourier transform of h. The authors expand h(u) in terms of Coiflets 
wavelets, 

+oo +oo +oo 

(5) fl(u) = 22 C m,k(pm,k( U ) + d hk^i,k( U )- 

k=—oo 3 = m k=— oo 

where m ,fc(<^) = 2 m / 2 0(2 m o> — A;), ipj :k (cu) = 2 J ' 2 ^(2 J 'o; — fc), being and ^ the scaling and 
wavelet functions respectively. 

The next step consists in inverting the expression (jSJ) by means of the Fourier inversion 
formula. Finally, considering the expression (pi}, the formulae of Laplace inversion become, 

fc=— oo ' 

/(x) = lim f m (x). 

m— >+oo 

One drawback of this approximation is that the wavelet approach involves an infinite 
product of complex series and the computation of the Fourier transform of some scaling 
functions. This can look intimidating for practical applications and may also take relatively 
long computational time. 

Based on operational matrices and Haar wavelets, the author in [1] presents a new method 
for performing numerical inversion of the Laplace transform where only matrix multiplications 
and ordinary algebraic operations are involved. However, the essential step in the method 
consists in expressing the Laplace transform in terms of -, which is impossible when we 
just know numerically the transform. Another drawback of this method is that the matrices 
become very large for larger scales. 

2. MULTIRESOLUTION ANALYSIS AND CARDINAL B-SPLINES 

A natural and convenient way to introduce wavelets is following the notion of multireso- 
lution analysis (MRA). Here we provide the basic definitions and properties regarding MRA 
and B-spline wavelets, for further information see [T3| 15]. 

Definition 1. A countable set {f n } of a Hilbert space is a Riesz basis if every element f of 
the space can be uniquely written as f = Y2 n c nfm o,nd there exist positive constants A and 
B such that, 

A\\f\\ 2 <Y,\ C n\ 2 <B\\f\\ 2 . 

n 

Definition 2. A function <p £ L 2 (IR) is called a scaling function, if the subspaces V m of 
L 2 (R), defined by, 

V m = clos L 2 {R) {<p mM : k e Z} , m e Z, 
where (f) m ,k — 2 m ^ 2 (j)(2 m x — k), satisfy the properties, 

(i) ■ ■ ■ C VLi C V Q C Vi C • • • . 

(ii) cl0S L i (\J m &Vm) =L 2 (R). 
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f^n me zKn = {o}. 

(iv) For each m, {4> m ,k ' k G Z} is a Riesz (or unconditional) basis ofV m . 
We also say that the scaling function 4> generates a multiresolution analysis {V m } of L 2 (M. 

The j th order cardinal 5-spline function, Nj(x), is defined recursively by a convolution: 
Nj(x) = / Nj-x(x - t)N (t)dt = / Nj_i(x-t)dt, j > 1, 



where, 

N (x) = X[o,i)(x) 

Alternatively, 



1 if x G [0, 1) 
otherwise. 



Nj(x) = %N j - 1 {x) + j + 1 X Nj-\(x - 1), j>l. 

We note that cardinal 5-spline functions are compactly supported, since the support of the 
jth orc [ er 5-gpline function Nj is [0, j + 1], and they have as the Fourier transform, 



IW 



In this paper we consider <fp = Nj as the scaling function which generates a MRA (see 
Figured]). Clearly, for j = we have the scaling function of the Haar wavelet system. We 
also remark that from the previous discussions, for every function f m G V m , there exists a 
unique sequence {c^^kez G / 2 (Z), such that, 

fcez 



Order 
Order 1 
Order 2 
Order 3 



Figure 1. Cardinal B-splines of orders j = 0, 1, 2, 3. 
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3. The wavelet approximation method 
Let us now consider a function / G L 2 (R) and its Fourier transform, whenever it exists: 

/+oo 
e- iwx f{x)dx. 
-oo 

Since / G L 2 (IR) we can expect that / decays to zero, so it can be well approximated in a 
finite interval [a, b] by, 



(6) f c (x) 



f(x) if x G [a, b], 
otherwise. 



Let us approximate f c (x) ~ fm,j( x ) f° r an x [ a i^]i where, 

0'+l).(2 m -l) , _ . 

(7) = + 

fc=0 ^ ' 

with convergence in L 2 -norm. Note that we are not considering the left and right boundary 
scaling functions (we refer the reader to Section 3 in [T2] for a detailed description of scaling 
functions on a bounded interval). 

The main idea behind the Wavelet Approximation method is to approximate / by /° • 
and then to compute the coefficients d m k by inverting the Fourier Transform. Proceeding 
this way, 



f(w)= / e™*f(x)dx~ / e- twx f^(x)dx 

J — oo J 
(i+l).(2--l) +oc 



-OO 



X) 



s ^(/r^(° +i) 'H)^) 

Introducing the change of variables y = (j + 1) • gives us, 

r_ C?+l)-(2™-l) „ +00 

fc=o - 7 - 00 

(j+1).(2™-1) . 

- 7 - , i e V 7 + 1 

J fc=0 VJ 

Finally, taking into account that (f> J mk (C) — ''4> : '(^Oe - *^ 5 "^ and performing the change of 

__ ! b—a -j. 

variables z = e 2m u+i) ( we have, 

(i+l)-(2" l -l) 

J ~k 
u m,k A 

k=0 

If we consider, 

(i+iH^-D 2f(j + l),- 2!! ^/(^±ili.log(,; 
Pm,j{z) = ) c^ mk z and Q m ,j(z) = - 

then, according to the previous formula, we have, 

(8) p m>j {z) ~ g mJ (^). 



~/2 m (j + l). , A ™6-a 2">o + i)a r , ' 
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Since P m j(z) is a polynomial, it is (in particular) analytic inside a disc of the complex plane 
{z G C : \z\ < r} for r > 0. We can obtain expressions for the coefficients c , mk by means of 
the Cauchy's integral formula. This is, 

■k-aS * = °.~.U + i)-(*"-i). 

where 7 denotes a circle of radius r, r > 0, about the origin. 

Considering now the change of variables z = re m , r > 0, gives us, 

(9) 

where A; = 0, (j + 1) • (2 m — 1), and and stand for the real and imaginary parts 
of 2, respectively. 

Note that if k 7^ then we can further expand the expression above by considering, 
(10) c4 ifc = [\(P m ,(re m ))cos(ku)du. 







On the other side, since ^ (i ■ log(z)) = ( j , we have, 

_ 2f (j + l),- 2 ^/ (gggi • log(s)) (log(,))^ 

(11) Q mj (z) - ______ , 

and it has a pole at z — 1. Finally, making use of (jSJ) and taking into account the former 
observation, we can exchange P m j by Q m ,j in @ and ([TP]) to obtain, respectively, 

(12) < * ^ J 2 \(Q m>j (rn)du, 
and, 

(13) - Jjf jT ^(Qm,i(re il1 )) cos(^)^, fc = 1, (j + 1) • (2 m - 1), 

where r 7^ 1 is a positive real number. 

In practice, both integrals in (1121) and ([13]) are computed by means of the Trapezoidal 
Rule, and we can define, 



I(k) = / U{Q m>j {re m )) cos(ku)du } 
Jo 

I{k] h) = ^ ( Q mJ (r) + (-l) fc Q m ,,(-r) + 2 ^ 3?(Q mJ (re^)) cos(£;/i s ) ) , 



3=1 

where h = and h s = sh for all s = 0, . . . , M. Proceeding this way we find, 

(14) 1 / M-l 

Q mj (r) + (-l) fc Q mj (- r ) + 2 £ M(Q m>j (re ih °)) cos(kh s 



Mr k 



where k = 1, (j + 1) ■ (2 m - 1). 

Let us summarize four sources of error in our procedure to compute the numerical Fourier 
transform inversion using cardinal B-splines wavelets. These are: 
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(A) Truncation of the integration range, 

£x(x) :=f(x)-f c (x), xeR\[a,b\. 

(B) The approximation error at scale m, 

£2(2) := f c (x) - f^ d (x), x G [a, b}. 

(C) The discretization error, which results when approximating the integral I(k) by I(k; h) 
using the trapezoidal rule. We can apply the formula for the error of the compound 
trapezoidal rule considering, 

O) = ^(QmA^ U )) cos(ku), £ 3 := I{k) - I(k; h), 

and assuming that q J mk G C 2 ([0,7r]). Then, 

(15) \^\ = ^p\(ikk(^yV J*e(o,7r). 

(D) The roundoff error. If we can calculate the sum in expression ffT^|) with a precision 
of 10~ v , then the roundoff error after multiplying by a factor is approximately 

jfpz -10 71 . Then, the roundoff error increases when r approaches to 0. 

3.1. Choice of the parameter r. As mentioned before, the choice of the parameter r may 
have influence in both, the discretization and the roundoff errors. In this section we present 
a detailed analysis of the errors listed before in order to determine the optimum r value. To 
do this, let us consider fi(x) = X[i/2,i)( x ) a step function defined in [0, 1], where x represents 
the indicator function, and its Fourier transform f±(w) 



Due to the shape of fi it seems that the best B-splines basis to perform the approximation is 
based in the Haar scaling functions (B-splines of order 0). In this particular case K(Q m j(re m )) 
and its derivatives up to order 2 can be computed relatively straightforwardly, so that we 
will be able to calculate the optimal value of the parameter r in order to minimize the 
discretization and the roundoff errors. We also demonstrate that the approximation error 
(type (B)) is in this case. To do so, we consider, 

q ( re ») _ r2 '" ( cos (2 mM ) + i sin(2 m w)) - r 2 " 1 " 1 (cos(2 m - 1 w) + % sin(2 m - 1 M)) 

2 m / 2 (r cos(m) +irsin(u) — 1) 

where u G [0,7r], and, 

r 2m+1 cos(2 m w - u) - r 2m cos(2 m -u) - r 2m ~ 1+1 cos(2 m ~ 1 w - u) + r 2 '"" 1 cos(2 m - 1 u) 
~ 2 m /2( r 2 _ 2rcos(n) + 1) ' 

Now, we must choose an appropriate r value to control both the discretization and the round- 
off errors. First of all, we consider the discretization error which can be estimated by means 
of expression (JT5|) . We note that q t ^ lk G C 2 ([0,7r]) since, 

< (r- l) 2 < r 2 - 2rcos(w) + 1 < (r + l) 2 , Mu G [0,7r],r ^ 1. 

So we have, 

= ^{Qm fl {re lu )) co S (ku) - m{Q mfl {re m )) sin(ku), 
au ' au 

-^q° m M = ^(Q m , {re iu )) cos(ku) - 2k^(Q mfi (re lu )) sm(ku) - k 2 U(Q m , (re iu )) cos(ku). 
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and, 

,2" ; -;-.. 

2 m / 2 



(16) \X(Q mfl (re-)\ < - + "J" - r 2 ^ + o(r 2 ^), 



where the last approximation holds for suitably small values of the parameter r. For sake 
of simplicity, we consider only the terms with smaller exponents in the parameter r for the 
expressions £9f?(Q mi o(W n )) and ^9f?(Q m , (>e m )). Then, 



(17) 
and, 
(18) 



2 (m- 2)/2 r 2- 1+A(r) 
^ fr _l)4 — ~ r + °( r < 



2 (3m-4)/2 r 2'"- 1 + 5 / r \ 

< -, 7TT — - ^ "' * + o(r 

(r — l) 8 



2?n — 1 , 2 m — 1 s 



where A(r) and B(r) are polynomials in r with degree greater than 2 m_1 , and the approxi- 
mations in (ITT)) and fjl8jl hold for suitably small values of the parameter r. Finally, taking 
into account expressions (TT5)) . f|T6|) . (ITTj) and ([TBI we have, 



^ " iW +fcV^ 1 )+o(r-- 1 ) = ^ (A; 2 + 2^ + l)r 2m - 1 +o(r 2m - 1 ). 

We note that j^l — > as r \, while the roundoff error is increasing in r as r approaches 
to zero. 

Then, the total error should be approximately minimized when the two estimates are equal. 
This leads to the equation, 

• 10-" = (k 2 + 2k + l) r 2m ~\ 

Mr k 12M 2 v ; 

After algebraic manipulation, we find, 

/ 12M-10~ r ' \ 

< 19 > ^= U* + tt + i) ) • * = ».-.s"-i- 

As mentioned before, the roundoff error arises when multiplying the sum in expression 
(fbil) by the pre-factor (Mr k )~ l . Let us consider M = 2 m , then the k of interest is k = 2 m — 1 
which is the greatest value that this parameter can take (small values of k do not cause 
roundoff errors). 

The left plot of Figure [2] represents the pre-factor for values of r > 0.9, while the right plot 
shows the pre-factor values for r > 0.999. We also display in Tabled] the pre-factor values 
(Mr k )~ l for different values of r and scales m = 8, m = 9 and m = 10. 

On the one hand, it is clear that we must concentrate in this second interval in order to 
get a reasonable roundoff error and on the other hand, the discretization error grows when r 
is close to 1. Later, in the numerical examples section, we will confirm the theory developed 
above for the step function. 

Since fi is compactly supported, we have E\ = 0. 

Note that in the case that j = 0, 

2 m -l , s +oo 2 l -l , s 

f (*) = £ ^ + £ £ H • 

k=0 ^ ' l=m k=0 ^ ' 
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r 


Scale 


m = 8 


m = 


9 


m — 


10 


0.9000 


1.8194 


■10 y 


4.7077 


1Q 2U 


6.3040 


10 43 


0.9100 


1.0869 


■ 10 8 


1.6618 


10 18 


7.7688 


10 38 


0.9200 


6.6968 


• 10 6 


6.2395 


10 15 


1.0833 


10 34 


0.9300 


4.2522 


• 10 5 


2.4885 


10 13 


1.7047 


10 29 


0.9400 


2.7807 


• 10 4 


1.0529 


10 11 


3.0193 


10 24 


0.9500 


1.8717 


• 10 3 


4.7203 


•TO 8 


6.0042 


10 19 


0.9600 


1.2960 


• 10 2 


2.2394 


• 10 6 


1.3373 


10 15 


0.9700 


9.2550 


■ 10° 


1.1230 


• 10 4 


3.3283 


10 10 


0.9800 


6.7470 • 


io- 1 


5.9457 


•10 1 


9.2347 


• 10 5 


0.9900 


5.0674 • 


io- 2 


3.3201 • 


IO" 1 


2.8503 


•TO 1 


0.9990 


5.0415 • 


10~ 3 


3.2566 • 


10~ 3 


2.7177- 


10~ 3 


0.9991 


4.9145 • 


10" 3 


3.0942 • 


IO" 3 


2.4532 ■ 


10~ 3 


0.9993 


4.6699 ■ 


10~ 3 


2.7934 ■ 


10~ 3 


1.9990- 


10~ 3 


0.9995 


4.4376 • 


10~ 3 


2.5219 • 


10~ 3 


1.6289- 


10~ 3 


0.9997 


4.2169- 


10~ 3 


2.2768 • 


10~ 3 


1.3274- 


10~ 3 


0.9999 


4.0072 • 


IO" 3 


2.0555 • 


IO" 3 


1.0818- 


IO" 3 



Table 1. Pre-factor (Mr k )~ l for M = 2 m , k = 2 m - 1 and scales m = 8, m = 9 and m = 10. 



where, fe }fc=o,...,2 m -i U \}i> m k=o,- ,2 l -i is the Haar basis (orthonormal) system in 

L 2 ([0,1]). Then, ' 

(20) 



l^2||iS([ a) 6]) — ||/ C frnW^da. 



b]) 



-too 2'-l , 
I—™ I — n n / 



l=m k=0 



+oo 2«-l 



since II* (f=f)IL 2(M) = (6- a) • KJ^^jj 
depends on the length of the interval [a, b] and the detail coefficients, 



2 l=m k=0 

b — a. Then, the approximation error 



(21) 



a i,k 



f c (x) ■ ipi tk (x)dx. 



Furthermore, since the detail coefficients ( 12T|) are zero, then we have also that £2 = and 
the approximation is exact at any scale level m. 
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4. The COS method 

For completeness, we present here the methodology developed in [7] for solving the inverse 
Fourier integraQ. The main idea is to reconstruct the whole integral from its Fourier- cosine 
series expansion extracting the series coefficients directly from the integrand. Fourier-cosine 
series expansions usually give an optimal approximation of functions with a finite support 
[3]. In fact, the cosine expansion of f(x) in x equals the Chebyshev series expansion of 
/(cos" 1 ^)) in t. 

For a function supported on [0,7r], the cosine expansion reads, 

A +oc 
fc=l 

with Ak = ~ f(6) cos(k9)d9. For functions supported on any other finite interval [a, b] e R, 
the Fourier-cosine series expansion can easily be obtained via a change of variables, 



It then reads 



with, 



x - a b - a n 

t) := - it, x = & + a. 

b — a 7r 



A x — a 

f( x ) = IT + Ak cos ( kn T )> 

2 ^-^ b — a 

k=l 



x — a s 



(23) A k = / f{x) cos(kiT- )dx. 

b — a J a b — a 

Since any real function has a cosine expansion when it is finitely supported, the derivation 
starts with a truncation of the infinite integration range in (|2"2"|) . Due to the conditions for 
the existence of a Fourier transform, the integrands in (1221) have to decay to zero at ±oo and 
we can truncate the integration range in a proper way without losing accuracy. 

Suppose [a, b] G K is chosen such that the truncated integral approximates the infinite 
counterpart very well, i.e., 

(24) CiW — f e iwx f{x)dx~ [ e lwx f{x)dx = £{w). 



Here, £i denotes a numerical approximation. 

Comparing equation (124"]) with the cosine series coefficients of f(x) on [a, b] in ( 123]) . we find 
that, 



where 9ft denotes the real part of the argument. It then follows from (1241) that A k — F k with, 



2 _ / , / kir 



K U e~ % —* 



; kciTT 



b — a 

We now replace A k by F k in the series expansion of f(x) on [a, b], i.e., 

A + °° — 
(25) h(x) = ^ + VF fc cos(A;7r^), 

2 ^-^ b — a 

k=i 



1 In order to maintain the notation used by the authors, here 
(22) £H = f e™*f{x)dx, 

represents the characteristic function, and hence the Fourier transform of a density function /(x) 
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and truncate the series summation such that, 

A N ~ l - 
(26) f 2 (x) = -f + J2F k cos(kn^-^). 

k=l ' 

The resulting error in f2{x) consists of two parts, a series truncation error from ( |25p to (1261) 
and an error originated from the approximation of A k by F k . An error analysis that takes 
these different approximations into account is presented in [7]. 

The COS method performs well when approximating smooth functions, but many terms 
are needed in case that the function or its first derivative has discontinuities along the domain 
of approximation. 

5. Numerical examples 

The aim here is to show the accuracy of B-splines to invert Fourier transforms of extremely 
peaked or discontinuous functions with finite support. We will consider the B-splines scaling 
functions of order j = 0, 1, 2. We denote by WAi-j the Wavelet Approximation method with 
B-splines of order % at scale j, which involves (2P(i + 1) — i) coefficients and COS-N the COS 
approximation method with N terms. 

To compute the coefficients of the Fourier inversion in expression (Tbil) . we must set the 
parameters. For this purpose we consider M = (j + 1) • 2 m , where j is the order of the 
B-spline considered and m is the scale parameter. Observe that if we take M = 2k instead 
of M = (j + 1) ■ 2 m in (TH1) this leads to the following expression, 

1 / 2fc_1 

s=l 

fc-1 

Q md (r) + (-l) fc Q m ,,(-r) + 2 ^(-l) s 3f?(Q mj (re^)) 



(27) 



2kr k 




2kr k 



for k — 1, + ■ (2 m — 1), so the computation time reduces for large scale approximations. 

Finally, we set r = 0.9995. Although we know that r > and r ^ 1, we must take into 
account the two types of errors (C) (discretization) and (D) (roundoff) listed previously, 
which may have influence on the computation of the coefficients. Due to the fact that the 
function ^R.(Q m j(re lu )) is in general intractable from an analytical point of view, we carried 
out intensive simulations in order to understand the influence of parameter r. We did these 
simulations considering the test functions fi and /^(x) = e - "'^, a > at different scale levels 
and used B-splines scaling functions of order j = 0, 1, 2. To show just an example, we have 
plotted the results for f 2 , a = 50 and B-splines of order 1 at scale m = 5 (Figure [3]) and m = 9 
(Figure H]). The colors represent the magnitude of the logarithm of the absolute errors. As 
we can observe, the absolute error remains constant for values \r — 1| < e. When \r — 1| > e, 
the error increases for high values of x (i.e. high values of the translation parameter k) and 
decreases for low values of x (i.e. low values of the translation parameter k). It is worth 
mentioning that for high scales, the error grows very rapidly when |r — 1| > e (empirical 
findings demonstrate that e ~ 0.05, for this reason we take approximately the midpoint of 
the interval avoiding r = 1), while this effect diminishes for shorter scales. In fact, as we will 
see later with the step function, the roundoff error almost disappears at very low scales (for 
instance m = 1) and the discretization error tends to zero when r tends to zero. These facts 
confirm the high impact of the roundoff error at high scale levels and only the impact of the 
discretization error at very low scales. 

For comparison, we gave in Section H] a brief overview of the COS method developed by 
Fang and Oosterlee ([?]) that is based on a Fourier-cosine series expansion and which usually 
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Figure 3. Logarithm of the absolute error for the approximation of the function /2 (a = 
50) with B-splines of order 1 at scale m = 5. 
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Figure 4. Logarithm of the absolute error for the approximation of the function /2 (a = 
50) with B-splines of order 1 at scale m = 9. 



gives optimal approximations of functions with finite support Q3J). The COS method is state 
of the art and has been applied to efficiently recover density functions from their Fourier 
transforms in order to solve important problems arising in Computational Finance. 

Let us consider the step function presented before. Observe that if we take m — 1, there 
is almost no roundoff error and the only remaining error is the discretization error £ 3 . If we 
assume that 77 = 16, then according to (DSD r 1>0 = 7.740368 • 10~ 17 , r M = 4.398968 ■ 1(T 9 . 
The left plot of Figure shows the approximation to the step function fx with the COS and 
the WA method with r = 0.9995. The right plot of Figure [5] represents the absolute error of 
the approximation with the COS and the WA method with r = 0.9995 and the optimum r 
computed with expression (TT9|) . Note that with just 2 coefficients the WA method is capable 
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to accurately approximate the step function, while the COS method suffers from the Gibbs 
phenomenon even if a lot of terms (we show up to 2048) are added to the COS expansion. 

Note that if we consider the approximation to the step function at scale m = 10, then 
the values tto^ computed for the parameter r are all of them in a neighborhood of 0.95 in 
accordance with the massive simulations performed before. 

The conclusion in that case is that the COS method performs poorly around the jump and 
exhibits the Gibbs phenomenon, while Haar wavelets are naturally capable to deal with this 
discontinuity. 




Figure 5. Zoom of the approximation (left) and absolute error of the approximation 
(right) to the function f\. 



5.1. Exponential function. Let us consider the exponential function, f%{x) = e a > 
and its Fourier transform f2(w) = a 2+ w 2 ■ 

Observe that this function becomes extremely peaked when a » 1. We will test the 
inversion with a = 50 and a = 500. 

We consider the interval [a, b] = [—1,1] and a = 50,500. We carry out the comparison 
between both methods using approximately the same number of terms in the approximation. 
It is worth mentioning that the computation of wavelet coefficients is more time consuming 
than the calculation of COS coefficients. Moreover, the COS method is easier to implement. 
Results are reported in Table [2] and plotted in Figure [6j We observe that linear B-splines are 
the most suitable basis functions to approximate highly peaked functions like the exponential 
we are considering. While adding many more terms in the COS expansion improves only a bit 
the approximation, when we consider B-splines of order 1 at higher scales the approximation 
error decays much more faster. The WA method with B-splines of order 1 performs better 
than the WA method with B-splines of orders and 2. 

In addition, we consider the following representative examples, 

5.2. B-spline basis function. We consider the function f%{x) = 2<p\ (x) in the interval 
M = [0,2]. 

We aim to recover the original function through the Wavelet Approximation inversion 
method. We apply also the COS method. Like in the previous example, the errors of type 
(A) and (B) do not have to be considered for this function, since is compactly supported 
and we carry out the approximation by means of B-splines of order 1. However, errors of 
type (C) and (D) remain. 

As before, we consider r = 0.9995. The recovered coefficients are c\ = 2, c\ 1 = —5.605010- 
10 -8 and c\ 2 = 4.645665 ■ 10~ 8 , showing high accuracy in the approximation. 
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WAO-6 


-2.928931 


-0.374819 


-2.523483 


-0.040001 
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-7.719564 
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-8.119854 


-0.083319 
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-12.704332 


-3.107169 


-9.290222 


-1.194033 
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COS-64 


-5.673413 


-0.526046 


-5.825399 


-0.057691 


COS-128 


-5.483920 


-0.805937 


-5.305269 


-0.120147 


COS-256 


-6.093322 


-1.102053 


-5.898013 


-0.244112 


COS-512 


-7.488787 


-1.402252 


-6.433619 


-0.450188 


COS-1024 


-7.454768 


-1.703285 


-6.500191 


-0.716606 



Table 2. Approximation errors to the function /2 in the interval [—1, 1]. 






COS-1024 
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Figure 6. Absolute error of the approximation to the function fa with a = 50 (top) and 
a = 500 (bottom). 



Figure [7] shows a zoom of the approximation with COS method to the function / 3 in a 
neighborhood of 0.5 with 64 and 128 terms, while the absolute error is plotted in the right 
part of the figure. As expected, the error grows up at the non-smooth part of the function, 
that is, near the points x = 0, x = 0.5 and x = 1. 

5.3. Combination of B-spline basis functions. We consider the function f^(x) = Ylt=o e 
1) in the interval [a, b] = [—1, 1], this is, a finite sum of basis B-splines of order one with coef- 
ficients which exhibit exponential decay. Again, the coefficients can be accurately recovered 
by the Wavelet Approximation method with B-splines of order one at scale 2. The maximum 
absolute error is, 
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Figure 7. Zoom of the approximation (left) and absolute error of the approximation 
(right) to the function fa. 



max \e k - c\ J = 2.681058 • lCT 8 . 

k=0+6 ' 

Left plot in Figure E] represents the approximation to the function with the COS method, 
while right plot shows the absolute error of the approximation in log scale. As before, the 
conflictive points for the approximation are the points of non differentiability. 

I 1 1 1 1 1B+000 r 1 1 1 1 




Figure 8. Approximation (left) and absolute error of the approximation (right) to the 
function f^. 



1 x 

5.4. Gaussian function. Finally, let us consider the Gaussian function, feyx) = ^^- e ^ 

-~- to 2 2 

and its Fourier transform f§(w) = e ~ . 

Since the function f$ is analytic, the COS method performs much more better than the 
Wavelet Approximation. We can see the results in Table [3j As expected, quadratic B-splines 
behave better than Haar or linear B-Splines, since the maximum error decreases when the 
order of the B-Spline increases. 

Left plot in Figure M shows the graph of the exponential function while the right plot 
represents the graph of the Gaussian function. It is worth observing the relation between the 
continuity of the function or its derivatives and the order of the B-Spline which fits better 
in the approximation. Haar wavelets for functions with a jump discontinuity, first order 
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B-splines for funtions with a jump discontinuity in the first derivative and finally, quadratic 
B- Splines for more regular functions. 



Method 


a 

min log | error 


= 0.1 

max log | error 


WAO-6 


-16.090813 


-0.428258 


WA1-5 


-11.687910 


-1.482131 


WA2-4 


-9.285967 


-2.450365 


COS-32 


-9.392864 


-5.390718 


COS-64 


-17.699992 


-7.530656 



Table 3. Approximation errors to the function /s in the interval [—1, 1]. 





Figure 9. The exponential function f2 with a = 50 (left) and the Gaussian funtion fa 
with a = 0.1 (right). 



6. Conclusions 

We have investigated a numerical procedure, the Wavelet Approximation method, to invert 
the Fourier transform of functions with finite support by means of the B-splines scaling 
functions. First of all, we truncate the function in a sufficiently wide interval and then 
we approximate it by a finite combination of B-splines wavelets up to order 2. Finally the 
function is recovered from its Fourier transform. 

We have tested and compared the accuracy of this method versus the COS method for a set 
of heterogeneous functions in terms of the continuity and smoothness. We have considered 
a continuous function with a jump discontinuity in its first derivative and a function with 
a jump discontinuity in its domain of definition. Additionally, combinations of B-splines 
of order one have been also considered as the functions to be recovered. With these last 
two examples we are able to assess the discretization and roundoff errors of the Wavelet 
Approximation method, since no errors of type (A) and (B) take place. Finally, we also have 
considered the Gaussian function. COS method performs better with infinitely differentiable 
functions, like the Gaussian, due to the fact that the coefficients of the expansion decay 
exponentially fast. On the contrary, WA method is more suitable for functions that exhibit 
peaks or discontinuities along its domain. For the function with a jump discontinuity, B- 
splines of order (Haar) are better, while B-splines of order 1 fit accurately the continuous 
function with a jump discontinuity in its first derivative. Furthermore, little improvement is 
achieved for these two last functions, when adding a lot of terms in the COS expansion. 
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Due to the fact that COS coefficients are easier to compute in terms of CPU time, hybrid 
methods involving COS and WA, or simply combinations of the WA method at different scales 
and/or with different orders of the B-splines, could be investigated in the future. B-splines are 
a semi-orthogonal system, so it would be also advisable to explore the orthogonal Daubechies 
or Battle-Lemarie scaling functions. Furthermore, it is well known that standardized B- 
splines tend to the Gaussian function when the order tends to infinity, so they could be very 
useful to recover Gaussian functions from its Fourier transform. 
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